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Abstract 

We present parton distribution functions which include a quanti- 
tative estimate of its uncertainties. The parton distribution functions 
are optimized with respect to deep inelastic proton data, expressing 
the uncertainties as a density measure over the functional space of 
parton distribution functions. This leads to a convenient method of 
propagating the parton distribution function uncertainties to new ob- 
servables, now expressing the uncertainty as a density in the prediction 
of the observable. New measurements can easily be included in the 
optimized sets as added weight functions to the density measure. Us- 
ing the optimized method nowhere in the analysis compromises have 
to be made with regard to the treatment of the uncertainties. 
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1 Introduction 



With the advent of new hadron colhder experiments the need to quan- 
titatively estimate parton distribution function (PDF) uncertainties is of 
paramount importance. One obvious example where the PDF uncertainty 
will become crucial is the ly-boson mass uncertainty. By no means is this the 
only observable for which the PDF uncertainty is of the utmost importance. 

The realization that hadron collider experiments had reached a level of 
accuracy where a quantitative approach towards the PDF uncertainties was 
needed came with the one-jet inclusive transverse energy distribution mea- 
surement of the CDF collaboration using the run la data with an integrated 
luminosity of 19.5 pb~^ |jl[]. Up to that point the qualitative "global fitting" 
approach had worked well and was used to demonstrate the success of the 
perturbative QCD framework. Be that as it may, the precision of the CDF 
one-jet inclusive measurement demonstrated that the qualitative approach 
had reached its limit in usefulness. The fact that prior to the measurement 
none of the PDF's could predict the high transverse momentum data whereas 
after the measurement the PDF's could be adjusted to accommodate the high 
transverse energy excess speaks for itself. The only way to continue is to 
use a quantitative approach which reflects the PDF uncertainties. 

The qualitative "global fitting" approach 0, |[ combines a large amount 
of experimental data to form an error weighted average of all input data, i.e. 
find the minimum solution. In the calculation of the the experimental 
statistical and systematic uncertainties are added in quadrature. Because 
experiments are combined with each other, even if no statistical adequate 
solution could be found, no conventional probabilistic interpretation can be 
reached (see for instance section II. D of ref. for an overview of the para- 
doxal problems encountered when attempting a probabilistic interpretation 
in the "global fitting" approach). 

It is crucial to develop a rigorous statistical approach such that the uncer- 
tainties have an objective statistical interpretation and in case of deviations 
the experiments or prior assumptions causing the conflict are traceable. The 
apparent way to proceed would be to redo the "global fitting" approach 
with the uncertainty analysis in mind. That is, treating the experimental 
error analysis more carefully by including correlation matrices reflecting the 
experimental uncertainty and by demanding the input experiments to be 
compatible. As an advantage this would still allow a calculation with a 
minimalization procedure. The final result is then quoted as a correlation 
matrix of the PDF parameters, i.e. a 2nd order Taylor expansion around the 
best solution. Yet, this method has problems. First of all the experimental 
systematic uncertainties are not gaussian despite the fact that often gaussian 
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approximations are published. Secondly, one can expect complicated corre- 
lations between the PDF parameters given the fact that the problem of PDF 
fitting is highly non-linear causing lower dimensional non-trivial regions in 
PDF parameter space with constant probability measure (i.e. constant x^)- 
Thirdly, in the end the PDF's are not the physical observables and prop- 
agation from the PDF parameters to the physical observable would again 
necessitate a linearization of the error propagation involving an estimate of 
the derivative of the observable with respect to the PDF parameters. In 
short, this method has three layers of subsequent linearizations of uncertain- 
ties which can lead to problems. 

The first paper using this method [P] uses the HI, ZEUS, BCDMS and 
NMC data. Subsequently, refs. |^ include additional experiments into the 
fits. Recently, a Lagrange multiplier method was outlined which does not 
need to make gaussian approximations 0. Also proposed in refs |TD| is 



a procedure to maintain the "global fitting" philosophy by include as many 
experiments as possible. To recover a probabilistic interpretation one has to 
multiply the experimental uncertainties with a common factor such that all 
experiments are consistent within a single PDF "global fit". 

Parallel to the previous methods, a different approach was pursued in 
ref. |jTl| to remove the last of the three gaussian approximations in the er- 



ror analysis. Using the gaussian parameter probability density of ref. an 
optimized Monte Carlo integration approach was applied to predict physical 
observables expressing its PDF uncertainty as a density probability measure 
in the predictions of the observable. Note that such an approach can be 
mathematically reformulated as a statistical inference method. In this paper 
we extend this method to eliminate all approximations in the error analysis. 
As a result any experimental error analysis can be implemented and the final 
resulting uncertainties for observables can assume any form. The interme- 
diate PDF parameter probability distribution can have any shape required 
by the experiments. The resulting optimized Monte Carlo approach is trans- 
parent, easy to use, modular and extendible. These properties are important 
for practical use. 

In section 2 we will explain our method of optimized Monte Carlo inte- 
gration. After that we are ready to optimize the PDF's towards sets of ex- 
periments. In section 3 we will consider the proton F2 results of HI, ZEUS, 
BCDMS and NMC together with an world averaged value of as{Mz)- The 
conclusions and outlook will be given in section 4. 
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2 The Method 



This section will describe the methodology we developed to obtain theoret- 
ical predictions of physical observables including PDF uncertainties. Before 
we can formulate our method we need to define carefully all aspects of com- 
paring experiments and theory. First of all, one has to define a theoretical 
framework which approximates the true nature value of observables. Sec- 
ondly, the experimental results have to be casted in a well defined object 
which will take the form of an experimental response function. After defin- 
ing the theoretical prior and the experimental response function we have the 
right language to formulate our method. 



2.1 The Model Priors 

Any observable has a true nature value xt which has to be calculated (i.e. 
approximated) by a theory model giving a value Xt- The PDF's only consist 
within the framework of the parton model and perturbative QCD . This 
implies approximations have been made to arrive at this calculable model. 
Hence the resulting PDF's will depend on the theory model defined in the 
prior. As long as the experimental uncertainties are larger than the deviations 
from the true nature value due to the theory model the approach is valid. 
The model uncertainties fall into three distinct classes. 

The first class is an agglomeration of non-perturbative effects. One has 
to make the approximation in which factorizable PDF's and calculable hard 



scattering matrix elements exist . This leads to the appearance of power 
suppressed momentum transfer terms. These higher twist terms can be either 
parametrized (see e.g. ref. |^) or assumed small enough to be neglectable by 
applying appropriate kinematic cuts. Moreover, the final state outgoing par- 
tons are assumed to have factorizable fragmentation functions. This again 
leads to power suppressed momentum transfer terms. Finally there is the in- 
teraction of the colliding hadron remnants with the hard scattering. For this 
effect no real physics models exist and usually the experiment is "corrected" 
for the "underlying event" fathering a large systematic uncertainty. For this 
paper we will neglect all such effects in the theory prior and assume they are 
neglectable with respect to the other uncertainties in play for the kinematic 
cuts applied. The underlying event subtraction is, as usual, considered an 
experimental problem. One final non-perturbative effect comes in play when 
experiments are considered which use other initial state hadrons than pro- 
tons (e.g. deuterium, iron, copper, etc.). In such a case one has to consider 
the effect of the multitude of nuclei. That is, what is the relation between the 
free proton PDF and the measured atomic densities. In this paper we will 



4 



not include heavy target data and therefore can ignore the phenomelogical 
shadow models which attempt to describe these nuclear effects. Note that 
the inclusion of non-perturbative models will involve the introduction of ad- 
ditional degrees of freedom (i.e. non-perturbative parameters). The effect is 
to de-emphasize certain kinematical regions associated with low momentum 
transfer scattering. Such non-perturbative models, if included in the theory 
prior, must be applied consistent to all predictions. 

The second class of uncertainties is related to the perturbative expansion 
of the hard scattering partonic cross sections. Because this is in principal 
calculable the uncertainties are more traceable. In the theory prior one has 
to determine the order in perturbative QCD to which all observables will be 
calculated. This also defines the order to which the PDF evolution has to 
be performed. A special case would be to define a resummation scheme of 
dominant initial state logarithms to extend the kinematic range of reliable 
predictions. If such a scheme is adopted it should be applied to all predic- 
tions. Also such a scheme could affect and modify the evolution equations. 
For this paper we choose the next-to-leading order approximations in the MS 
renormalization/factorization scheme for all observables. For the appropriate 
evolution of the PDF's the program QCDNUM was used [|13| which accuracy is 



more than adequate for current phenomenology |14]. Note a theoretical un- 



certainty could be assigned to the fixed order hard scattering matrix element 
calculation which reflects an estimate of the deviation from the prediction 
to an even higher order calculation of the observable. This can take many 
forms and will be highly subjective. However, one can define such a prior 
and incorporate it easily in the uncertainty analysis. In this paper we restrict 
ourselves to only looking at the renormalization/factorization scale depen- 
dence of the next-to-leading order theory predictions. To do this we allow a 
floating renormalization/factorization scale in the optimalization procedure. 
Given enough orders in the perturbative expansion virtually no scale de- 
pendence would remain. However, at next-to-leading order there will be a 
scale dependence on the observables and therefore the PDF's. The extend 
to which the scale can vary independent of the PDF probability measure 
indicates the need to increase the perturbative order of the matrix element 
calculation. Given enough experimental accuracy the fixed order calculation 
will fail. This will be reflected through the scale in a strong preference for a 
specific value of this scale. 

The final class of model priors is a set of additional requirements. Some of 
them can be physics necessities, e.g. the PDF charge and momentum conser- 
vation sum rules. Others can be more speculative, for instance assumptions 
about small and large parton fraction behavior or introducing PDF moment 
constraints from lattice QCD (see e.g. WW^)- Also, one could use a previ- 
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ous PDF optimalization as a prior and include additional experiments in the 
optimized PDF sets. Note that this requires that the prior optimalization re- 
sults are consistent with the new experiments, i.e. that the inclusion leads to 
a refinement of the relevant region contributing to the functional integration 
in the functional space V{J^) of PDF's JF. The most important constraint 
in this class are smoothness constraints on the actual PDF's which must be 
introduced to regulate unconstrained fiuctuations due to the discrete nature 
of experimental results. The most restrained introduction of some of these 
requirement is to define a specific parameterization for the functional form 
of the PDF's depending on a fixed, finite number of parameters. Note how- 
ever that this still requires an assumption on the initial probability density 
distribution of those parameters. For this paper we choose one such parame- 
terization as is detailed in the next section. All in all, the constraints have to 
be quantified and assembled as a prior probability density function measure, 
Ppriori^), over the PDF functional space V{J^). 

2.2 The Experimental Response Function 

Given the theory prior we now can make well defined approximations Xt of 
the true nature value xt for any relevant observable. In order to interact 
with the experimental results we have to define the experimental response 
function 

-Pea;p('^e I "^('^t ) ) Pexpi^-^ e\-^ o) ; (1) 

which is a probability density estimating the likelyhood of measuring Xe given 
the theory prediction xt and detector response V. The detector response 
is the deformation infiicted on the physics signal by the detector resulting 
in the measured signal Xq = V{xt). The systematic uncertainties are now 
defined as a probability density measure, Pgysi'D), over the functional space 
V{V) of all potential detector responses for the observable and quantifies 
the understanding of the detector. Note that in principle the experimental 
response function has to be formulated before the actual measurement of the 
particular observable. Also, theory predictions can originate from any source 
and are not tied to a specific model in any way. 

Often the detector response is specific to a particular measurement and 
not correlated with any other experiment under consideration. In that case 
we can integrate over the detector response functional space 

PexpiXe\Xt) = j dV Pexp{Xe\V{xt)) X Psys{V) , (2) 

thereby simplifying the analysis. For instance, if the detector response prob- 
ability function is parameterizable in a multi-gaussian the above integration 
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would simplify to the usual formulation with a correlation matrix encap- 
sulating the systematic uncertainties. 

Note that if several measurements with correlated systematic uncertain- 
ties are included one has to integrate the systematic uncertainties over the 
group of correlated measurements. A more flexible approach is to use the 
optimized Monte Carlo approach (explained in the next subsection) not only 
for the PDF's but in conjunction with the detector response T>. This way, 
each PDF has a specific detector response D which makes propagation of 
systematic uncertainties to other observables trivial. 

The experimental response function Pgxp together with the detector re- 
sponse function Pgyg contains all information we can extract from the mea- 
surement. Together with the actual measured values of the observable, Xm, 
they form a permanent and well defined record of the measurement. 



2.3 The Optimized Monte Carlo Approach 

We now can formulate our method. To make a prediction for observable O 
which includes the PDF uncertainty based on prior assumptions and certain 
sets of measurements, we perform the integration over the PDF and detec- 
tor response functional space to obtain the probability density function of 
observing a value Xe for observable O 

where P^^^ (xglxf (JF)) is the experimental response function for observable O 
using PDF set T. Note we have integrated out the systematic uncertainties 
for the observable as explained in eq. |^ and thereby assuming its systematic 
uncertainties are independent of the other experiments involved in the PDF 
determination. The prior probability function Pprior is defined in subsection 
2.1 and quantifies all prior PDF assumptions and regulators. Finally, the 
probability function of the input experiments is defined as 

(4) 

with the individual experimental response function for each experiment 

Pe.p k{J',V) = Pe,.p k {xl^^\V (xf ^(^))) , (5) 

where the actual measurement is substituted in the experimental re- 
sponse function. While formally this defines the PDF uncertainty for any 
observable we want, in practice a numerical integration over the PDF's and 
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possible systematic uncertainties has to be implemented. The only feasible 
method in this situation is a Monte Carlo approach. That is, we approxi- 
mate the integral by randomly picking a sufficient large set of PDF's and 
possible detector response functions such that 

O ^ ^ inut o 

and 

Pexp (*^ i) ~ Pexp iy T^i) X Pexp 2(-^ j; T^i) X ■ ■ ■ X Pexp ri(-^ T^i) X PgysiT^i) ■ 

(7) 

While this can be numerically implemented, the Monte Carlo integration 
efficiency defined as the ratio of the average of the combined probability 
function and the maximum of the combined probability function 

e(P(^,P)) ^ (^EP(-^., A)) / (maxP(^„I?,)) , (8) 

with 

P{J^,V) = Pe^p i{J^,V)xPe^p 2(J^,P)X- • -xPe^p n{.^,V)xP,y,{V)xPp,.ior{.^) , 

(9) 

is for all practical purposes zero. This because given the virtually infinite 
functional space of choices of JF and T) the likelihood of choosing at random 
a set which gives a non-neglectable contribution to the functional integral is 
zero. 

The standard method to tackle this problem is to optimize the Monte 
Carlo integration with respect to the combined probability function, i.e. "un- 
weighting" . Looking at eq. |^ this procedure simply comes down to changing 
the integration measure by redefining the PDF's JF to JF' together with the 
detector response function V to V such that all the probability densities 
are absorbed in the integration measure. After such a transformation eqs. ^ 
and H simplify to 




Translating this in the Monte Carlo integration approach now gives an effi- 
ciency of one, i.e. the Monte Carlo sets are "unweighted" and eqs. |^ and ^ 
transform in the trivial formula 

1 ^ 
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where the probabihty measure is quantified in the density of func- 
tional sets. 

Note that we could also have optimized with respect to a subset of the 
experiments, e.g. experiments 2 through n. We would then find 




which tells us how to include an additional measurement, not included in the 
optimalization procedure. However, care has to be taken with the efficiency. 
The effective number of PDF's used to estimate the observable is no longer 
N but iVg// = e{Pexp i) X N . If the efficiency is too low the Monte Carlo esti- 
mate uncertainty on the observable will become large and dominate over the 
PDF uncertainties. This can happen if the newly included experiment has a 
superior uncertainty analysis compared to the previous included experiments 
in which case one has to increase the initial number of sets N . Or, alterna- 
tively, the newly included experiment disagrees with the previous included 
experiments in which case the situation cannot be resolved by increasing A^. 
Coercing the experiments together plainly leads to an incorrect result. 

2.4 Comparison Methodology 

As explained in section 2.3 the predictions are based on a discreet set of 
PDF's and will build up a probability density in the space of the observ- 
able as is expressed in eq. |n|. Note that the probability function density 
function of measuring Xe for observable C, Pp^jr{xe), depends on the experi- 
mental response function of the measurement. Each discreet prediction will 
be "smeared" into a continuous probability function and averaged over all 
PDF's in the set. This means that the particular form of the experimen- 
tal response function dictates the minimum number of PDF's needed for a 
satisfactory result. A drawback is that each experiment measuring the ob- 
servable will have a different P^f- In more theoretical studies one often wants 
a prediction independent of any particular experiment. In that case we have 
to substitute for the experimental response function an idealized detector 
model. The most straightforward model is the "perfect" detector, i.e. 

1 ^ 

^pd/(a^e)^^E^(^e-X,(^.)) (13) 

where the sum runs over the optimized PDF's in the set. Such a result 
leads to a traditional scatter plot representation of P^j- By introducing a 
averaging parameter A one can get to a more continues result. One could 
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for instance use a gaussian with a width A for the response function. How- 
ever, traditionally one chooses two theta functions to get a histogrammed 
representation of P^j. 

1^1 1 
PpZfixe) ^ ^ E ©(^e - 2 A - X e(x,(^.) - - -A) . (14) 

1=1 

The drawback of such procedures is that resulting approximation of P^j- 
depends on the averaging procedure. 

While the above procedures will give us an approximation for the prob- 
ability density function of the observable it does not represent a likelyhood 
or confidence level. Using P^j(xe) it is pretty straightforward to construct 
the confidence level. However some care has to be taken given the potential 
non-gaussian nature of the experimental response function. First we define 
the parameter independent log-likelyhood 

L\X,) = -2 X l0g(P5/Xe)/n^,_,(Xe)) , (15) 

where 

Pexpmaxi^e) = I^^-XP (Xel^t) , (16) 

is the maximum obtainable probability given an observed value Xe- For a 
gaussian response function the log-likelyhood defined here resorts back to the 
usual definition. The definition of log likelyhood as a ratio of probability 
densities has the property that it is independent under reparameterizations 
of the observable. We now can define the confidence level 

CL<^(x,) = jdx QiL\x) - L\x,)) X P^/x) (17) 

which is a proper probability. Given the measured value Xm, the confidence 
level CL[xm) quantifies the probability that a repeat of the experiment will 
have a worse agreement with the theory prediction including the PDF un- 
certainties. 



3 The Optimized Sets 

The previous section defined the theory prior with the exception of the actual 
PDF parameterization. In subsection 3.1 the parameterization used for this 
paper will be discussed. Once that is done we perform in subsection 3.2 
the actual optimalization using various combinations of proton and as 
measurements. 
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3.1 The Parameterization Choice 



Given this is our first venture into the refractory subject of PDF determina- 
tion we choose the well established MRST-parameterization |jl6[| as our first 



guidance. This choice implicitly incorporates into the theory prior the accu- 
mulated knowledge of many years of PDF studies which can be considered 
a positive feature. However, because we have to adhere to the experimen- 
tal systematic uncertainties and statistical interpretation of the results this 
parameterization can turn out to be far too restrictive. Such restrictive 
parameterizations might result in at least an underestimate of the PDF un- 
certainties. Or, more seriously, in discrepancies between the theory and the 
data stemming not from physics but from the parameterization choice. See 



for instance ref. [[T^] where the gluon PDF parameterization effects on the 
one jet inclusive transverse energy distribution are discussed. 

For completeness, the explicit MRST-parameterization at the scale Qq = 
1 GeV is given by 

xg = AgX-^^{l-xy'^{l + eg^/^ + -fgx) 
xS = 2x{d + u + s) = Asx-^%l-xy^{l + es^/x + -fsx) 

xA = x{d-u) = Aax~^^{1- xy-''^^{l + -fAX + 5ax'^) . (18) 

The normalization coefficients Au, Ad and Ag are determined by the charge 
and momentum conservation sum rules. Both the charm and bottom quark 
PDF's are generated through perturbative evolution from mass threshold. 
For the charm quark the threshold mass is set to 1.5 GeV, while the bottom 
quark threshold mass is chosen to be 4.5 GeV. 

To conclude, using the MRST-parameterization the functional PDF inte- 
gration of subsection 2.3 is reduced to an integration over 21 parameters 
(not counting the strong coupling constant as and squared renormaliza- 
tion/factorization scale yu). The prior probability distribution of the pa- 
rameters is chosen to be uniform. Next, we have to optimize the parameters 
according to the combined probability density function. 



3.2 The rr and Data Optimized Sets 

Given the choice to abstain for the moment from any non-perturbative mod- 
eling in the theory prior, we have to limit ourselves to proton target data. 
The study of heavy nuclei data by introducing a shadowing model is a sub- 
sequent step in the development of the PDF's and will be taken in another 
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Experiment 


Measurement 


usable points 


used a:-range 


error 


analysis 


BCDMS 


pP 


344 


7.0 X 10-2-0.75 




gaussian 


HI 


pP 
^2 


188 


5.0 X 10-^-0.32 


half 


gaussian 


ZEUS 


pP 
^2 


187 


6.3 X 10-^-0.51 




gaussian 


NMC 


pP 
^2 


127 


8.0 X 10-3-0.50 




gaussian 


E665 


pP 
^2 


53 


3.7 X 10-3-0.39 




gaussian 


LEP 


as 


1 


N.A. 




gaussian 



Table 1: A summary of relevant properties of experiments considered in this 
paper. 



paper. The five deep inelastic F2 experiments we have selected so far are 
BCDMS H], NMC HI [H, ZEUS ^ and E665 [||. For the as mea- 



surement we take the value from ref [23|: as{Mz) = 0.119 ± 0.003. While 



this value is in fact a world average of sorts, we will indicate it in the fits 
as the "LEP" experiment. All the experiments are summarized in table 1 
together with the relevant properties. 

With the exception of HI, all other experiments quote a gaussianized 
analysis of the systematic uncertainties (i.e. a correlation matrix of the ex- 
perimental errors is either quoted or can be calculated using the published 
results). The HI F2 measurement has 5 sources of systematic uncertainties 
quoted as ±lcr deviations where for some of the sources the 0"+ is differ- 
ent from the (T_. This means that for HI we use the technique described 
in the previous section and numerically integrate over the systematic un- 
certainties by optimizing the PDF parameters together with the 5 detector 
model parameters of HI using as the Pgysi'D) detector response function 5 
half-gaussians depending on the detector model parameters. For the other 
4 experiments the correlation matrix is used, i.e. the systematic uncertain- 
ties are integrated out analytically. This implies all potential correlations 
between the different experiments are ignored. 

Because no higher twist models are in the theory prior we apply a cut on 
the momentum transfer and parton fractions in the F.[ data as suggested in 
ref. 0: g2 > 2 GeV^ and = x {1/x - 1) > 10 GeV^. 

Apart from optimizing the PDF integration over the 21+1 PDF param- 
eters and, for HI, the 5 detector model parameters, we also include the 
squared renormalization/factorization scale as an optimalization parameter. 
As explained in subsection 2.1 the reason is twofold. First of all varying the 
scale amplifies the PDF uncertainty stemming from the fixed order pertur- 
bative expansion approximation. Secondly, the resulting optimized factor- 
ization/renormalization scale will indicate to which degree the theory model 
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is applicable. A very narrow distribution indicating a failure of the model, 
a wide distribution indicating the model is appropriate. Note that we have 
chosen the initial squared scale probability distribution uniform. 



To obtain an optimized set we use a Metropolis Algorithm com- 
bined with a simulated annealing procedure^ [|25|. We will discuss some of 
its potential shortcomings as they can affect the results. The first potential 
problem is the possibility of equivalent disconnected maxima. While the an- 
nealing should guard against finding local maxima and nudge the Metropolis 
walker to the global maximum it can happen that more than one maximum 
exists of roughly equal probability which are not connected by a likely path. 
Given the finite number of steps the Metropolis walker takes this can cause 
problems as the generated optimized PDF's are concentrated on one of the 
maxima neglecting the others. A trail-and-error method of choosing different 
starting points for the Metropolis walker can be employed. However for a 
finite number of steps this does not guarantee the absence of a secondary 
maximum. The second problem arises from the choice of step size for the 
Metropolis walker. One wants to keep the step size small relative to the 
width of the probability function in order to maintain a good efficiency for 
generating optimized sets. However, this leads to correlations between PDF's 
subsequently generated. A common sense approach is to pick a small ran- 
domly ordered subset out of the generated optimized PDF parameter points. 
As that might be, for complicated topologies, especially subspace regions of 
constant probability, a large number of PDF's can be required to sample the 
entire region important to the functional integral. 

Note that these potential problems in generating the optimized sets are 
numerical in nature and not a shortcoming of the employed method. Both 
problems disappear in the limit of an infinite number of steps. Any finite 
number of Metropolis steps leaves the potential for problems as described 
above. For now we use the FERMILAB PC-farm to generate optimized 
sets of 100,000 PDF's (i.e. PDF parameters). Out of the 100,000 sets smaller 
subsets of 100, 1,000 and 10,000 PDF's are constructed by random selection 
from the larger list. For most practical applications 100 optimized PDF's are 
already enough. For more complicated analysis additional optimized PDF's 
might be necessary. 

The first optimalization to consider are with respect to the individual 
experiments. Given a set optimized to one of the experiments we can 
calculate the confidence level that the prediction using this set describe the 
other four experiments. The confidence level, as is detailed in section 2.4, 



^The temperature in the annealing procedure is identical to the tolerance parameter of 
ref. 1^, A value larger than 1 means the experimental uncertainties are amplified. 
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HI 


BCDMS 


E665 


ZEUS 


NMC 


LEP 


Hl-MRST set 




67% 


21% 


0.5% 


<0.1% 


31% 


BCDMS-MRST set 


85% 




23% 


1.5% 


<0.1% 


0.5% 


E665-MRST set 


30% 


82% 




1.6% 


1.0% 


99% 


ZEUS-MRST set 


22% 


<0.1% 


5.0% 




<0.1% 


24% 


NMC-MRST set 


<0.1% 


28% 


1.5% 


<0.1% 




3.2% 



Table 2: The confidence level of each experiment given the different sets. 
The name of the set is composed of all included experiments and the PDF 
parameterization choice. 



is defined as the probability a repeat of the experiment would result in a 
worse agreement than the current agreement. Given this analysis we can 
conclude which experiments are compatible and can be considered in com- 
bined optimized sets. The results are listed in table 2. Note that because 
different experiments will have different experimental response functions and 
hence different PDF uncertainties the confidence level of experiment A given 
experiment B is different from the confidence level of experiment B given 
experiment A. A dramatic point in case can be seen in table 2 between the 
HI and ZEUS experiment. From the results exposed in table 2 we conclude 
that both ZEUS and NMC cannot be combined with any other experiment 
as they produce their own distinct PDF sets. Both NMC and BCDMS pre- 
fer a much lower value of as than the LEP result. We can combine HI, 
BCDMS and E665 as each of these experiments predicts the other two with 
high confidence. 

Given the confidence level results we have constructed, apart from the 5 
sets using the individual F2 data results, 7 additional PDF sets with different 
combinations of HI, BCDMS, E665 and LEP results. Note that we combine 
the as measurement with BCDMS even though the confidence level tells us 
it is not compatible. While inconsistent we do this to get PDF sets with a 
somewhat more stabilized and acceptable as- Alternatively, we could have 
chosen to fix the value of as to a predetermined value in the physics prior. 

With the twelve optimized sets we are now ready for comparisons with 
hadron collider data. Given this is the first kind of such an endeavor we can 
expect surprises which will guide us to the next step in the further develop- 
ment of PDF's with uncertainties. In subsequent papers we will explore the 
phenomelogical implications for many hadron collider observables. 
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Figure 1: The Hl-MRST set parameter distributions (histograms) together 
with the gaussian approximation (red curve) and the central MRST fit resuh 
||T6| (blue bar). In addition the log-likelyhood AL^ distribution (histogram) 
and the probability distribution approximation (red curve) is shown. Note 
that the last 5 parameters are the detector model parameters of HI 
together with the prior probability distribution for this parameter (green 
curve). The green curve for the as result is the LEP value 
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3.3 Some PDF set Results 



Given the factorization scheme dependence required for the calculation of 
the next-to-leading order matrix elements, the individual PDF's do not con- 
tain much physics. Moreover the correlations between different PDF's can 
be expected to be significant given the fact that at leading order the F2 
experiments only determine the squared electric charge weighted sum of the 
flavor PDF's . At next-to-leading order one can hope at best to distinguish 
this charged sum of quark PDF's from the gluon PDF. This renders a metic- 
ulous study of a specific flavor PDF prodigal. Yet, a cursory look might 
be helpful to understand the phenomelogical implications for hadron collider 
observables. 

First of all we can look at the 21 PDF parameters of an optimized set 
together with as, the squared renormalization/factorization scale and the 
log-likelyhood distribution. In fig. 1 we take as an example the Hl-MRST 
optimized set which also includes the five detector modeling parameters. 
Prior to the optimalization each detector modeling parameters is given by 
a gaussian random parameter with an average of zero and a width of one. 
To obtain the fractional detector deformation of the signal a positive 
value of the random variable is multiplied by 0"+ and a negative value by a-. 
Examining fig. 1 we make several general observations which apply also to 
the other optimized set results: 

• Many of the parameters deviate significantly from gaussian behavior, 
letting the linearized error propagation method of refs. |^, |[ go 

astray. 

• The central MRST fit is not that far off for each individual pa- 
rameter. The differences are attributed to the different treatment of 
the systematic errors and the MRST fit is averaged over many more 
experiments. 

• The optimized detector response parameters are notably close to the 
HI central value, though some of the distributions are skewed. 

Also of interest are parameters not directly related to the PDF's but im- 
portant for the quality of the optimized set. These are the log-likelyhood 
distribution, the squared renormalization scale distribution and the as dis- 
tribution for each of the twelve sets. 

The value of as can be used to infer the reliability of the experimental 
error analysis. The confidence level intervals of the as distribution for all 



16 





a5-distribution 


AL^-distribution 


set 


CL interval 


mm(AL2) 


(AL2) 


io-2(AL2) 


NDP/DOF 


ZEUS-MRST 


(0.113, (0.114,(0.115), 0.116),0.117) 


466.1 


481.2 


42.9 


187/23 


NMC-MRST 


(0.098, (0.102,(0.108), 0.112),0.117) 


185.4 


196.4 


10.1 


127/23 


Hl-MRST 


(0.108, (0.112,(0.115), 0.117),0.119) 


166.6 


175.9 


7.8 


188/28 


Hl+LEP-MRST 


(0.114,(0.116,(0.118), 0.119),0.121) 


167.3 


176.0 


8.3 


189/28 


BCDMS-MRST 


(0.104,(0.106,(0.108), 0.110),0.112) 


317.2 


328.1 


12.1 


344/23 


BCDMS+LEP-MRST 


(0.112, (0.113,(0.116),0.117),0.119) 


325.0 


335.8 


15.7 


345/23 


E665-MRST 


(0.106, (0.112,(0.116), 0.127),0.133) 


57.9 


65.5 


4.9 


53/23 


E665+LEP-MRST 


(0.114,(0.117,(0.120), 0.123),0.126) 


59.1 


66.5 


6.0 


54/23 


Hl+BCDMS-MRST 


(0.109, (0.110,(0.112),0.114),0.115) 


510.9 


525.8 


11.4 


532/28 


Hl+BCDMS+LEP-MRST 


(0.110,(0.111, (0.112),0.114),0.115) 


511.5 


521.8 


10.0 


533/28 


H1+BCDMS+E665-MRST 


(0.109, (0.111, (0.112), 0.114),0.115) 


580.3 


596.2 


12.3 


585/28 


H1+BCDMS+E665+LEP-MRST 


(0.110,(0.112,(0.113),0.114),0.115) 


579.7 


592.3 


10.4 


586/28 



Table 3: The relevant properties of the as and AL^ distributions for the 
optimized sets. The confidence level intervals are for a CL of 4.55%, 31.73% 
and 100%. The bin width used to calculate the confidence level intervals is 
0.005 using the 1,000 PDF's 



twelve PDF sets are listed in table 3. The notation used, {a, {b, {c),d),e), 
is the following: the 100% confidence level interval is given by c (i.e. the 
maximum value); the 31.73% confidence level interval is given by the interval 
(&, d) (i.e. the "1-sigma" interval); the 4.55% confidence level interval is given 
by the interval (a, e) (i.e. the "2-sigma" interval). In order to calculate the 



confidence level we use eq. with given by eq. [13: 

1 N 

CL{as) = T7 E (Ppdfi^^s (^.)) - Ppdfias)) • (19) 
i=i 

To evaluate the Ppdf functions in the step function we have to measure the 
density of as predictions. To do this we choose the histogram approxima- 
tion of eq. |14| with a width of 0.005. The number of PDF's involved in the 
evaluation, A^, is chosen to be 1,000. We can see the inclusion of as into the 
HI and BCDMS optimalization pulls as towards the LEP value. For HI it 
more than halves the uncertainty on as, while for BCDMS the uncertainty 
is not affected much. This is understandable because the BCDMS result 
and the LEP result excluded each other up to a confidence level of 0.5%. 
Hence, coercing them into a combined optimalization has shifted the opti- 
malization region and not, as is the case for the HI optimalization, refined 
the optimalization region. For the sets combining several measurements 
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we see that the inclusion of the LEP data point hardly affects the already 
well established value of as- This can be understood by the large number of 
datapoints included in the combined experiments, rendering the uncertainty 
on as significantly smaller than the uncertainty on the LEP value of as- 

From the log-likelyhood we can deduct something about the quality of 
the optimalization. The log-likelyhood distribution is for all twelve sets close 
to the expected distribution. The relevant parameters are given in table 
3. The minimum log-likelyhood, min(AL^), found is on average equal to the 
number of data points (NDP) minus the degrees of freedom (DOF). How- 
ever, due to the high degree of correlation between the PDF parameters the 
effective degrees of freedom (i.e. the independent degrees of freedom taking 
the correlations into account) is often much lower. The effective degrees of 
freedom are on average given by half the squared variance, |cr^(AL^). The 
difference between the average value of the log-likelyhood, (AL^), and the 
minimum value is for a distributions equal to the effective number of in- 
dependent parameters. From table 3 we can see that the ZEUS set and to a 
lesser extend the NMC set are problematic. This does not come as a surprise 
given the confidence level results of table 2. 



Using the renormalization/factorization scale dependence one can deduce 
the applicability of the theory model. All 12 distributions are shown in fig. 2. 
The troubled ZEUS and NMC sets are also exposed here by the narrow range 
of optimized squared scales signaling a distress between the experimental 
result and the theory model. Yet, some of the other distributions are also 
worrisome indicating possibly a too restrictive parameterization choice. Note 
that while in the prior we assumed an uniform probability distribution for 
the squared scale, we restricted the range to be larger than half the squared 
momentum transfer as we have a kinematic cut of > 2 GeV^ and a 
parameterization scale of Qo = 1 GeV^. 

Finally in fig. 3 we compare the optimized sets with the MRS99 0] and 
CTEQ5M 0] distributions. This comparison forecasts the deviations we will 
see in the hadron collider observables between the global fitter predictions 
and the optimized PDF sets. As can be seen the squared charged sum of 
the flavored PDF's agrees pretty well with MRS99 and CTEQ5M up to the 
point where there is data. Some subtle differences are present in some of 
the sets. The big difference however is the gluon distribution which is much 
lower above a parton fraction of 0.1 and higher between 0.01 and 0.1. These 
differences are consistent for all optimized PDF sets using the F2 data and 
will be reflected in the phenomenology at hadron colliders. 
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Figure 2: The squared renormalization/f actor ization scale distribution of the 
optimized sets. The horizontal axis is the squared ratio of the scale and the 
momentum transfer. 
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It is interesting to note that the differences in the gluon PDF are especially 
large when compared with CTEQ5M. This gives us a first hint at the origin of 
differences with the optimized sets. The relevant difference between MRS99 
and CTEQ5M is the choice of experiment to constrain the gluon PDF in the 
large x region. The MRS99 global fit makes the traditional choice of using 
the WA70 prompt photon data. Because the prompt photon data suffers 
from large theoretical uncertainties the CTEQ5M global fit instead chose the 
one jet inclusive transverse jet energy data of CDF and DO |2^. Yet the 
deep inelastic proton data favors the MRS99 over the CTEQ5M global fit 
result for large x gluons. This is a rather troubling conclusion as it suggests 
a very significant deviation between the deep inelastic proton data and high 
transverse energy jet data at the TEVATRON. 
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Figure 3: The ratio of the optimized PDF sets (blue) and CTEQ5M (green) 
over MRS99 at a scale of 10 GeV. The dashed red curves are the two MRS99 
sets with as{Mz) values of 0.1125 and 0.1225. The left column is the charged 
summed combination of PDF's , while the right column ratio of the gluon 
PDF . 
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4 Conclusions and Outlook 



In this paper we introduced a flexible method to incorporate PDF uncertain- 
ties into phenomelogical predictions of collider observables. By representing 
the uncertainties as a density, i.e. an ensemble of PDF's with which one pre- 
dicts the observable, the practical use is both convenient and flexible. The 
sets used in this paper can be obtained from the website pdf .fnal.gov. 

The restriction to a minimum set of experiments to determine the opti- 
mized PDF's is deliberate as many experiments do not seem to be consistent 
with a single optimized PDF set. One of the next steps in developing the 
optimized PDF's is to consider deuterium data and the related nuclear ef- 
fects in the optimalization. While for deuterium the nuclear effects are small 
and manageable within the context of PDF optimalization, one can expect 
serious complications if one considers heavier targets. It might be better to 
use the optimized PDF's to predict these measurements in order to obtain 
better understanding of the nuclear effects. 

In the current optimized sets the most severe prior assumption is due to 
the restrictive MRST parameterization. One can argue the seemingly dis- 
crepancies between proton F2 measurements and high energy jet production 
at CDF and DO is due to the restrictive parameterization. It is therefore cru- 
cial to device more general schemes for parameterizing the PDF functionals. 
One way is to use complete sets of functions, in such a way that the higher 
order functions are associated with smaller scale fluctuation. The drawback 
of such methods is that changes in a parameter affect the PDF for all values 
of the parton fractions x. This induces strong correlations between the pa- 
rameters and more importantly accurate data in localized regions of x will 
affect the uncertainty estimates for all regions of x. A better method would 
be splines or other methods of generating functions using the value of the 
PDF at fixed grid points Xi. In such a case the parameters (i.e. the value 
at the grid points) are more localized. Such methods are under development 
and the subject of future publications. 

Using the current optimized PDF sets a series of studies on hadron col- 
lider observables will be published. These studies range from luminosity 
measurements using H^-boson and Z-boson cross sections to jet production. 
Given the upcoming run II at Fermilab and the LHC program at Cern fur- 
ther development of the optimized PDF's is more than warranted as it will 
have a large impact on the phenomenology at these hadron colliders. 
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